FVCOM_R_Subsetting

Using Bounding Polygons to Trim FVCOM Mesh Data

Published

March 25, 2024

Accessing FVCOM in R

This approach will lean on the work of Ben Tupper at Bigelow Labs, and their package fvcom. Many thanks to Ben and their team for their work on this.

Read more about their package here: https://github.com/BigelowLab/fvcom

Accessing NECOFS GOM4

As explained on the repo documentation for Ben’s {fvcom} package, Data are served via OpeNDAP on a THREDDS Server. NECOFS has its own THREDDS Directory which can be used to browse what is available.

Working with the Nodes & Elements

We will use this Jan, 2019 data to get information about the Nodes we wish to gather data from.

[1] "2019-01-01 00:00:00 UTC" "2019-01-01 01:01:52 UTC"
[3] "2019-01-01 01:58:07 UTC" "2019-01-01 03:00:00 UTC"
[5] "2019-01-01 04:01:52 UTC" "2019-01-01 04:58:07 UTC"

Working with the Mesh

Finding the Nodes We Care About:

We can do ourselves a favor by identifying which nodes we’re interested in before pulling variables. There are a couple different approaches to try:

Use the nodes as points, use st_intersect to identify which are within the areas we are studying.
2. Use the mesh itself, to preserve the triangular geometries. Then use the {fvcom} functions to pull data

Areas We are Interested in

VTS Survey, Nearshore Lobster Habitat
GOM & GB Ecological Production Units
Northeast Shelf

Mesh Containment Within an Area

Will test this first with the offshore areas. But plan is to take the mesh and check for overlap/containment within the above polygons

Coordinate Reference System:
  User input: EPSG:6318 
  wkt:
GEOGCRS["NAD83(2011)",
    DATUM["NAD83 (National Spatial Reference System 2011)",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Puerto Rico - onshore and offshore. United States (USA) onshore and offshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore and offshore."],
        BBOX[14.92,167.65,74.71,-63.88]],
    ID["EPSG",6318]]

Use New Mesh to Extract Vars:

Extracting Vector Elements

Interpolate Currents to Regular Grid

Next Steps - GOM3

To go back further than 2016 we need to use the gom3 hindcast records. These have a different file location on the THREDDS directory, and we’ll need to confirm the grid structure and variables are consistent across these different model runs.

Verify that the IDS for nodes and elements are the same between NECOFS and FVCOMGom3.

If we want to go back further in time then NECOFS won’t yet be online and we’ll need to use data from the FVCOM GOM3 hindcast. Best-case scenario these mesh id’s match so we can append the time-periods together

Can we verify they match? Good.

What about the time frequency, what is the pace of the hindcast data records? Looks hourly, cool thats probably too much.

[1] 31.26466

Save Node/Element ID’s for regions of interest

Now that we know that the node and element ID’s are consistent across GOM3 & NECOFS, we can use the approach above to save the mesh information for areas within the regions we are looking at:

These can be used if-needed to limit the data being transferred over the OPeNDAP connection and save us data transfer/processing time. Would also be useful to have just as a means to match points to a standalone mesh geometry.

Coordinate Reference System:
  User input: EPSG:6318 
  wkt:
GEOGCRS["NAD83(2011)",
    DATUM["NAD83 (National Spatial Reference System 2011)",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Puerto Rico - onshore and offshore. United States (USA) onshore and offshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore and offshore."],
        BBOX[14.92,167.65,74.71,-63.88]],
    ID["EPSG",6318]]

Download data for an extended period of time.

We don’t need hourly, or even daily records. But we will need to make sure we’re aggregating in time correctly and not just grabbing snapshots.

So where we sit now. We know the exact ID’s for the nodes and centroids that make complete mesh elements that fall within our area of interest. We aslo know the variables we are interested in, and possibly the time period we want to cover and the temporal resolution we’re interested in (months).

$names
 [1] "name"          "len"           "unlim"         "group_index"  
 [5] "group_id"      "id"            "dimvarid"      "units"        
 [9] "vals"          "create_dimvar"

$class
[1] "ncdim4"

$names
 [1] "name"          "len"           "unlim"         "group_index"  
 [5] "group_id"      "id"            "dimvarid"      "units"        
 [9] "vals"          "create_dimvar"

$class
[1] "ncdim4"

$names
 [1] "nprocs"        "partition"     "x"             "y"            
 [5] "lon"           "lat"           "xc"            "yc"           
 [9] "lonc"          "latc"          "h"             "nv"           
[13] "nbe"           "ntsn"          "nbsn"          "ntve"         
[17] "nbve"          "a1u"           "a2u"           "aw0"          
[21] "awx"           "awy"           "art2"          "art1"         
[25] "iint"          "Itime"         "Itime2"        "Times"        
[29] "zeta"          "file_date"     "u"             "v"            
[33] "omega"         "ww"            "ua"            "va"           
[37] "temp"          "salinity"      "km"            "kh"           
[41] "kq"            "q2"            "q2l"           "l"            
[45] "short_wave"    "net_heat_flux" "uwind_stress"  "vwind_stress" 
[49] "fvcom_mesh"

$hasatt
[1] FALSE

$value
[1] 0

$hasatt
[1] FALSE

$value
[1] 0

[1] 48451    46

[1] -1.000001  0.000000

[1] 99137

start count 
    1    99

[1] 6

--- title: "FVCOM_R_Subsetting" description: | Using Bounding Polygons to Trim FVCOM Mesh Data date: "Updated on: `r Sys.Date()`" format: html: code-fold: true code-tools: true df-print: kable self-contained: true editor: visual execute: echo: false warning: false message: false fig.height: 6 fig.width: 8 fig.align: "center" comment: "" --- ## Accessing FVCOM in R This approach will lean on the work of Ben Tupper at Bigelow Labs, and their package fvcom. Many thanks to Ben and their team for their work on this. Read more about their package here: https://github.com/BigelowLab/fvcom ```{r} #| eval: false # To install fvcom #devtools::install_github("BigelowLab/fvcom") ``` ```{r} library(sf) library(fvcom) library(ncdf4) library(tidyverse) library(gmRi) ``` ### Accessing NECOFS GOM4 As explained on the repo documentation for [Ben's {fvcom} package](https://github.com/BigelowLab/fvcom?tab=readme-ov-file), Data are served via OpeNDAP on a THREDDS Server. NECOFS has its own [THREDDS Directory](http://www.smast.umassd.edu:8080/thredds/catalog/models/fvcom/NECOFS/catalog.html) which can be used to browse what is available. ```{r} # We can Access links from the THREDDS directory as if they are NETCDF files uri_base <- "http://www.smast.umassd.edu:8080/thredds/dodsC/models/fvcom/NECOFS/Archive/NECOFS_GOM/2019" uri <- file.path(uri_base, "gom4_201901.nc") x <- nc_open(uri) ``` ### Working with the Nodes & Elements We will use this Jan, 2019 data to get information about the Nodes we wish to gather data from. ```{r} # Nodes and Elements x_nodes <- fvcom_nodes(x, what = "lonlat") x_elems <- fvcom_elems(x, what = "lonlat") # Time information, don't need it for this but good to know: head(fvcom_time(x, internal = TRUE), n = 6) # Can extract variables using node or element sst <- get_node_var(x, var = "temp") ``` ### Working with the Mesh ```{r} # Can also get the mesh itself as a simple feature collection mesh <- get_mesh_geometry(x, what = 'lonlat') # And we can use the mesh to request variables with it plot(sf::st_geometry(mesh), border = scales::alpha("gray", 0.6)) ``` ### Finding the Nodes We Care About: We can do ourselves a favor by identifying which nodes we're interested in before pulling variables. There are a couple different approaches to try: 1. Use the nodes as points, use st_intersect to identify which are within the areas we are studying. \ **2. Use the mesh itself, to preserve the triangular geometries. Then use the {fvcom} functions to pull data** # Areas We are Interested in 1. VTS Survey, Nearshore Lobster Habitat 2. GOM & GB Ecological Production Units 3. Northeast Shelf ```{r} # Read things in proj_path <- cs_path("mills", "Projects/Lobster ECOL") vts_poly <- read_sf(str_c(proj_path, "Ecological Data/Spatial_Boundaries/VTSsurvey_nearshore_area.geojson")) res_shapes <- cs_path("res", "Shapefiles") epu_path <- str_c(res_shapes, "EPU/") gom_poly <- read_sf(str_c(epu_path, "individual_epus/GOM.geojson")) gb_poly <- read_sf(str_c(epu_path, "individual_epus/GB.geojson")) shelf_poly <- read_sf(str_c(epu_path, "EPU_extended.shp")) st_crs(shelf_poly) <- st_crs(gom_poly) ``` ## Mesh Containment Within an Area Will test this first with the offshore areas. But plan is to take the mesh and check for overlap/containment within the above polygons ```{r} # Prep CRS st_crs(mesh) gom_poly <- st_transform(gom_poly, st_crs(mesh)) # Flag the locations that are within the domain gom_mesh <- mesh %>% st_join(gom_poly, join = st_within) %>% drop_na(EPU) %>% dplyr::select(-c(EPU, Shape_Leng, Shape_Area, source)) ggplot() + geom_sf(data = gom_mesh, aes(color = "FVCOM Mesh")) + geom_sf(data = gom_poly, aes(color = "Gulf of Maine EPU"), fill = "transparent") + labs("Mesh Triangles Within Gulf of Maine") ``` ## Use New Mesh to Extract Vars: ```{r} # These are the core variables at the nodes node_var_list <- c('temp', 'salinity') # Then Can we Grab what we want? gom_vars <- get_mesh( x, # Dataset lazyloaded with ncdf4 from THREDDS y = 1, # integer, or indices for siglay or siglev (depth indices) vars = c(node_var_list), # Variables we want mesh = gom_mesh, # Mesh to get them for time = 24 # 24th Time interval ) # Plot Them ggplot(gom_vars) + geom_sf(aes(fill = temp), color = "white") + scale_fill_distiller(palette = "RdBu") + map_theme() ``` ## Extracting Vector Elements ```{r} # Information at the centroids elem_var_list <- c('u', 'v') #Northward and Eastward Water Velocity # Use the unique elem ID's to get data at the points we care about gom_vects <- get_elem_var( x, # Dataset lazyloaded with ncdf4 from THREDDS y = 1, # integer, or indices for siglay or siglev (depth indices) var = elem_var_list, elem = unique(gom_mesh$elem), time = 1 ) # Or locate the element centroid locations and use them to get the relevant lon/lat info for this stuff gom_elems <- fvcom_elems(x, what = "lonlat", index = unique(gom_mesh$elem)) %>% left_join(gom_vects) %>% rename(clon = lon, clat = lat) # Join them back to the mesh gom_vars <- left_join(gom_vars, gom_elems) # Plot the Vector information as vectors on a map # Looks messy b/c its not a regular grid library(metR) ggplot(gom_vars) + geom_sf(aes(fill = temp), color = "transparent") + geom_vector( aes(x = clon, y = clat, dx = u, dy = v), arrow.angle = 30, arrow.type = "open", arrow.length = .25, pivot = 0, preserve.dir = TRUE, direction = "ccw", color = "gray20", alpha = 0.25)+ scale_fill_distiller(palette = "RdBu") + map_theme() ``` ### Interpolate Currents to Regular Grid ```{r} gom_elems_regular <- raster::stack( sapply( #c(node_var_list, elem_var_list), elem_var_list, function(f) { fvcom::rasterize(gom_vars, field = f) }, simplify = FALSE)) # Use Rastervis to make vectorplot library(rasterVis) rasterVis::vectorplot( gom_elems_regular, isField = 'dXY', col.arrows = "white", main = 'Surface Currents') ``` # Next Steps - GOM3 To go back further than 2016 we need to use the gom3 hindcast records. These have a different file location on the THREDDS directory, and we'll need to confirm the grid structure and variables are consistent across these different model runs. #### Verify that the IDS for nodes and elements are the same between NECOFS and FVCOMGom3. \ If we want to go back further in time then NECOFS won't yet be online and we'll need to use data from the FVCOM GOM3 hindcast. Best-case scenario these mesh id's match so we can append the time-periods together ```{r} # The hindcast has a differnent URL in the THREDDS Catalog # source: http://www.smast.umassd.edu:8080/thredds/hindcasts.html?dataset=fvcom/hindcasts/30yr_gom3 gom3_url = "http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3" gom3_x <- nc_open(gom3_url) # Get the mesh itself as a simple feature collection to compare against the other gom3_mesh <- get_mesh_geometry(x, what = 'lonlat') # And we can use the mesh to request variables with it plot(sf::st_geometry(gom3_mesh), border = scales::alpha("gray", 0.6), main = "GOM3 30-year HindCast Mesh") ``` Can we verify they match? Good. ```{r} # Check they match, from an eyeball check it looks like it, yay, even with the geomtries # st_drop_geometry(mesh) == st_drop_geometry(gom3_mesh) ``` What about the time frequency, what is the pace of the hindcast data records? Looks hourly, cool thats probably too much. ```{r} # 30 years, in daily records fvcom_time(gom3_x) %>% as.Date() %>% length() / (365 * 30) ``` #### Save Node/Element ID's for regions of interest Now that we know that the node and element ID's are consistent across GOM3 & NECOFS, we can use the approach above to save the mesh information for areas within the regions we are looking at: These can be used if-needed to limit the data being transferred over the OPeNDAP connection and save us data transfer/processing time. Would also be useful to have just as a means to match points to a standalone mesh geometry. ```{r} # Prep CRS st_crs(mesh) vts_poly <- st_transform(vts_poly, st_crs(mesh)) gom_poly <- st_transform(gom_poly, st_crs(mesh)) gb_poly <- st_transform(gb_poly, st_crs(mesh)) # Flag the locations that are within the domains mesh_trim <- function(mesh, domain){ # Will remove these columns after st_join names_drop <- names(domain)[which(names(domain)!="geometry")] new_mesh <- mesh %>% st_join(domain, join = st_within) %>% drop_na() %>% dplyr::select(-all_of(names_drop)) return(new_mesh) } # Run all Four vts_mesh <- mesh_trim(mesh = mesh, domain = vts_poly) gom_mesh <- mesh_trim(mesh = mesh, domain = gom_poly) gb_mesh <- mesh_trim(mesh = mesh, domain = gb_poly) #shelf_mesh <- mesh_trim(mesh = mesh, domain = shelf_poly) # Check them ggplot() + geom_sf(data = vts_mesh, aes(color = "mesh"))+ geom_sf(data = vts_poly, aes(color = "bounds"), fill = NA) # Save them fvcom_out <- str_c(proj_path, "FVCOM_support/") # # Do we want to save the geometry with? Guess so... # st_write(vts_mesh, str_c(fvcom_out, "VTSSurvey_FVCOM_nodes.geojson"), overwrite = T) # st_write(gom_mesh, str_c(fvcom_out, "GOM_EPU_FVCOM_nodes.geojson"), overwrite = T) # st_write(gb_mesh, str_c(fvcom_out, "GB_EPU_FVCOM_nodes.geojson"), overwrite = T) ``` ```{r} # How do they look when loaded vts_check <- read_sf(str_c(fvcom_out, "VTSSurvey_FVCOM_nodes.geojson")) # Lookin Goood! ggplot() + geom_sf(data = vts_check, aes(color = "mesh"))+ geom_sf(data = vts_poly, aes(color = "bounds"), fill = NA) ``` #### Download data for an extended period of time. We don't need hourly, or even daily records. But we will need to make sure we're aggregating in time correctly and not just grabbing snapshots. So where we sit now. We know the exact ID's for the nodes and centroids that make complete mesh elements that fall within our area of interest. We aslo know the variables we are interested in, and possibly the time period we want to cover and the temporal resolution we're interested in (months). ```{r} # VTS mesh is somehow smallest so lets use that one: vts_mesh <- read_sf(str_c(fvcom_out, "VTSSurvey_FVCOM_nodes.geojson")) # Lets use the hindcast as our access point: gom3_url = "http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3" gom3_x <- nc_open(gom3_url) # Do some digging to see if we can throw integers in for depth... attributes(gom3_x$dim$siglev) attributes(gom3_x$dim$siglay) attributes(gom3_x$var) # What depth layers do we want depths <- ncvar_get(gom3_x, "siglev") ncatt_get(gom3_x, "siglev", "units") ncatt_get(gom3_x, "siglay", "units") # There are 46? depth (siglev), and they exist at all the centroids, # but their spacing is not consistent (stretches for deep areas, compresses for shallow) dim(depths) range(depths) # Doesn't really match up length(unique(mesh$elem)) # Time index navigation: fvcom_nav( x = gom3_x, vals = c(1, 100), #vals = c(as.POSIXct("1978-01-01 00:00:00 UTC"), as.POSIXct("1978-06-01 00:00:00 UTC")), what = "time") # Manual lookup which(fvcom_time(gom3_x) == "1978-01-01 00:00:00 UTC") ```